home *** CD-ROM | disk | FTP | other *** search
/ Trading on the Edge / Trading On The Edge - CD-ROM Toolkit (Wayzata Technology)(2031)(1994).bin / mac / Mac_Files / Software Utilities / NN PreProcessing / FINNLDS.C < prev    next >
C/C++ Source or Header  |  1992-09-30  |  27KB  |  749 lines
  1. /* 15:14 09-20-92  (finnlds.c)  Non-linear Dynamical System Analysis */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6.  
  7. /************************************************************************
  8.  * Copyright(C) 1992 High-Tech Communications.                          *
  9.  * 103 Buckskin Court, Sewickley, PA 15143                              *
  10.  *                                                                      *
  11.  * Written by Casimir C. Klimasauskas                                   *
  12.  *                                                                      *
  13.  * All rights reserved.  No part of this program may be reproduced,     *
  14.  * stored in a retrieval system, or tramsmitted, in any form or by any  *
  15.  * means, electronic, mechanical, photocopying, recording or otherwise  *
  16.  * without the prior written permission of the copyright owner,         *
  17.  * High-Tech Communications.                                            *
  18.  *                                                                      *
  19.  * These programs are supplied on an "as-is" basis with no warranties   *
  20.  * of fitness or operability, either express or implied.                *
  21.  *                                                                      *
  22.  ************************************************************************
  23.  */
  24.  
  25. #define ASof(x) (sizeof(x)/sizeof(x[0]))
  26. #define    EMIN (1.e-15)
  27.  
  28. /* The routines in this module were translated by or created by
  29.  *   Casimir C. Klimasauskas.  Due to a lack of reference data and
  30.  *   the numerical intensity of the computations, no guarantee can
  31.  *   be made that the procedures are correct.
  32.  */
  33.  
  34. /************************************************************************
  35.  *                                    *
  36.  *            Interface Definintion                *
  37.  *                                    *
  38.  ************************************************************************
  39.     These routines have been interfaced to the utilities developed and
  40.     described in the previous issue of Advanced Technology for Developers.
  41.     From the main menu:
  42.  
  43.       hurst                - enter the "hurst" analysis section
  44.       series = series.label,start,end - enter a new series
  45.       startwindowsize = 2        - initial window size
  46.       deltawindowsize = 2        - amount to increment window size
  47.       iterations = 20        - number of iterations (max 100)
  48.       derivative = no        - use derivative of series
  49.       printfile = untitled.prn,type    - print output file name
  50.                       type: tab, comma, text
  51.       type [file]            - type a file to the screen
  52.       go                 - execute the command
  53.       quit                - return to main menu system
  54.  
  55.       correlationdim             - enter correlation dim
  56.       series = series.label,start,end - enter a new series
  57.       embedding = 5            - embedding dimension
  58.       tau = 2            - lag for phase space
  59.       {      radius = 0.1            - initial radius    }
  60.       {      delta = 0.1            - increment for radius    }
  61.       iterations = 20        - number of iteraitons (max 100)
  62.       printfile = untitled.prn,type    - print output file
  63.                       type: tab, comma, text
  64.       type [file]            - type a file to the screen
  65.       go [startdim,enddim]        - run with embedding =start..end
  66.       quit
  67.  
  68.       lyapunov                 - enter lyapunov section
  69.       series = series.label,start,end - enter a new series
  70.       embedding = 5            - embedding dimension
  71.       tau = 2            - lag for phase space
  72.       minscale = 0.1        - minimum radius
  73.       maxscale = 1.5        - maximum radius
  74.       evolve = 3            - number of time steps per evolution
  75.       lagminimum = 2        - minimum lag between points
  76.       printfile = untitled.prn,type    - print output file
  77.                       type: tab, comma, text
  78.       type [file]            - type a file to the screen
  79.       go [startdim,enddim]        - run with embedding =start..end
  80.       quit
  81.  */
  82.  
  83. /************************************************************************
  84.  *                                    *
  85.  *    Utilities Requried in Computation of Various Parameters        *
  86.  *                                    *
  87.  ************************************************************************
  88.  */
  89.  
  90. /* --- AverageD() --- Compute the average of a series of data points */
  91.  
  92. static double AverageD( vFP, n )    /* compute the average of a series */
  93. float        *vFP;        /* pointer to floating vector */
  94. int         n;        /* number of items in list */
  95. {   double     v=0.;        /* work accumulator */
  96.     int         i=n;        /* work index */
  97.     while( --i >= 0 ) v += *vFP++;
  98.     return( n>0? (v/((double)n)):0. );
  99. }
  100.  
  101. /* --- StdDevD() --- Compute the standard deviation of data points */
  102.  
  103. static double StdDevD( vFP, n )        /* compute the standard deviation */
  104. float        *vFP;        /* pointer to vector */
  105. int         n;        /* # of items in series */
  106. {   double     XBarD=0.;    /* average */
  107.     double     DevD=0.;    /* standard deviation */
  108.     double     vD;        /* work value */
  109.     int         i;        /* work counter */
  110.     float    *FP;        /* work float pointer */
  111.     if ( n > 1 ) {
  112.     /* NOTE: use this method to compute standard deviation to
  113.      * reduce problems with round-off errors when the average is
  114.      * very different from zero.
  115.      */
  116.     for( FP=vFP, i=n; --i >= 0; ) XBarD += *FP++;
  117.     XBarD /= (double)n;
  118.     for( FP=vFP, i=n; --i >= 0; ) {
  119.         vD = (XBarD - *FP++);
  120.         DevD += vD * vD;
  121.     }
  122.     DevD = sqrt( DevD / n );
  123.     return( DevD );
  124.     }
  125.     return( 0. );
  126. }
  127.  
  128. /* --- RegressI() --- Perform linear regression on two series --- */
  129.  
  130. static int RegressI( mDP, bDP, xFP, yFP, n )    /* compute regression line */
  131. double        *mDP;        /* slope */
  132. double        *bDP;        /* intercept:  y = mx + b */
  133. float        *xFP;        /* x-vector */
  134. float        *yFP;        /* y-vector */
  135. int         n;        /* length of vector */
  136. {   double     XBarD=0.;    /* x-average */
  137.     double     YBarD=0.;    /* y-average */
  138.     double     NumD=0.;    /* numerator */
  139.     double     DenD=0.;    /* denominator */
  140.     double     v;        /* work value */
  141.     double     m, b;        /* slope, intercept */
  142.     register float *aFP, *bFP;    /* work pointers */
  143.     register int i;        /* work counter */
  144.  
  145.     /* NOTE: this method is used to compute the regression equations,
  146.      * because it is much more stable than the numerically "quicker"
  147.      * methods if the mean for X & Y differ significantly from zero.
  148.      */
  149.  
  150.     for( aFP=xFP,bFP=yFP,i=n; --i >= 0; ) {
  151.     XBarD += *aFP++; YBarD += *bFP++;
  152.     }
  153.  
  154.     /* --- handle too few points for valid regression --- */
  155.  
  156.     if ( n < 2 ) { m = 0.0; b = YBarD; goto Done; }
  157.  
  158.     XBarD /= (double)n;
  159.     YBarD /= (double)n;
  160.  
  161.     for( aFP=xFP,bFP=yFP,i=n; --i >= 0; ) {
  162.     v = *aFP++ - XBarD;
  163.     NumD += v * (*bFP++ - YBarD);
  164.     DenD += v * v;
  165.     }
  166.     m = NumD / (DenD > EMIN? DenD:EMIN );
  167.     b = YBarD - m * XBarD;
  168.  
  169. Done:
  170.     if ( mDP ) *mDP = m;
  171.     if ( bDP ) *bDP = b;
  172.     return( 0 );
  173. }
  174.  
  175. /************************************************************************
  176.  *                                    *
  177.  *    HurstI() - Hurst Coefficient Computations            *
  178.  *                                    *
  179.  ************************************************************************
  180.     NOTE: the return vectors are always of length "IterN" long.  They
  181.     must be allocated at least this length.
  182.  
  183.     hurstFP[0] = Hurst coefficient over all windows
  184.     hurstFP[1] = Hurst coefficient over windows 0,1,2
  185.     hurstFP[2] = Hurst coefficient over windows 1,2,3
  186.     ...
  187.  */
  188.  
  189. #define    HREG    5    /* number of items for regression line */
  190.  
  191. int HurstI( hurstFP, nFP, rsFP, lognFP, logrsFP,
  192.         xFP, nxI, StartN, DeltaN, IterN, DiffFI )
  193. float        *hurstFP;    /* (r) hurst coefficients */
  194. float        *nFP;        /* (r) number of windows */
  195. float        *rsFP;        /* (r) R/S raw */
  196. float        *lognFP;    /* (r) log (number of windows) (x) */
  197. float        *logrsFP;    /* (r) log (R/S) (y) */
  198. float        *xFP;        /* (i) time-series as a vector */
  199. int         nxI;        /* (i) # of items in time series */
  200. int         StartN;    /* (i) starting window size */
  201. int         DeltaN;    /* (i) amount to increase window size */
  202. int         IterN;        /* (i) # of iterations */
  203. int         DiffFI;    /* (i) =0, use raw data; =1, use diffs */
  204. {
  205.     float    *wxFP=0;    /* work vector */
  206.     float    *rswFP=0;    /* work RS vector (rscolumn+6) */
  207.     double     MD, SD, RSD;    /* M=Average; S=Std Dev; RSD=R/S */
  208.     double     MnD, MxD;    /* min & max */
  209.     double     SumXtD;    /* total Xt (running) */
  210.     double     SlopeD;    /* slope of the regression line */
  211.     double     InterceptD;    /* intercept of the regression line */
  212.     float     vF;        /* floating work value */
  213.     int         CounterI;    /* current iteration */
  214.     int         WindowSizeI;    /* current window size */
  215.     int         NoOfWindowsI;    /* number of windows */
  216.     int         WindI;        /* window number */
  217.     float    *FP;        /* work float pointer */
  218.     int         wxI;        /* work index */
  219.     int         rcI=0;        /* return code */
  220.  
  221.     /* --- clear out results area --- */
  222.  
  223.     for( wxI = 0; wxI < IterN; wxI++ ) {
  224.     nFP[wxI]     = 0.;
  225.     rsFP[wxI]    = 0.;
  226.     lognFP[wxI]  = 0.;
  227.     logrsFP[wxI] = 0.;
  228.     hurstFP[wxI] = 0.;
  229.     }
  230.  
  231.     /* --- allocate dynamic memory --- */
  232.  
  233.     wxFP  = (float *)malloc( ((unsigned)nxI+2) * sizeof(float) );
  234.     rswFP = (float *)malloc( ((unsigned)nxI+2) * sizeof(float) );
  235.     if ( wxFP == (float *)0 || rswFP == (float *)0 ) {
  236.     rcI = -1;
  237.     goto Done;
  238.     }
  239.  
  240.     /* --- copy over the series or the difference --- */
  241.  
  242.     if ( DiffFI ) nxI--;        /* one less item if diffs */
  243.     for( wxI = 0; wxI < nxI; wxI++ ) {
  244.     if ( DiffFI ) {
  245.         if ( xFP[wxI] < EMIN || xFP[wxI+1] < EMIN ) vF = 1.;
  246.         else vF = log(xFP[wxI+1]) / log(xFP[wxI]);
  247.     } else vF = xFP[wxI];
  248.     wxFP[wxI] = vF;
  249.     }
  250.  
  251.     /* --- go through each iteration --- */
  252.  
  253.     WindowSizeI = StartN;        /* starting window size */
  254.     for( CounterI = 0; CounterI < IterN; CounterI++ ) {
  255.  
  256.     /* --- Compute the number of windows --- */
  257.  
  258.     NoOfWindowsI = nxI / WindowSizeI;    /* # of windows */
  259.     if ( NoOfWindowsI < 1 ) break;        /* nothing more to do */
  260.  
  261.     /* --- process each window --- */
  262.  
  263.     for( WindI = 0; WindI < NoOfWindowsI; WindI++ ) {
  264.         FP = &wxFP[ WindI * WindowSizeI ];    /* start of this window */
  265.         MD = AverageD( FP, WindowSizeI );    /* Average */
  266.         SD = StdDevD(  FP, WindowSizeI );    /* Standard Deviation */
  267.         SumXtD = 0.;
  268.         for( wxI = 0; wxI < WindowSizeI; wxI++ ) {
  269.         SumXtD += (FP[wxI] - MD);
  270.         if ( wxI == 0 || SumXtD < MnD ) MnD = SumXtD;
  271.         if ( wxI == 0 || SumXtD > MxD ) MxD = SumXtD;
  272.         }
  273.         RSD = (MxD - MnD) / (SD < EMIN? EMIN:SD);
  274.         rswFP[WindI] = RSD;        /* R/S for this window */
  275.     }
  276.  
  277.     /* --- Compute RS --- */
  278.  
  279.     RSD = AverageD( rswFP, NoOfWindowsI );
  280.     rsFP[   CounterI] = RSD;
  281.     nFP[    CounterI] = WindowSizeI;
  282.     logrsFP[CounterI] = log( RSD );
  283.     lognFP[ CounterI] = log( (double)WindowSizeI );
  284.  
  285.     /* --- Compute the next window size --- */
  286.  
  287.     if ( DeltaN == 0 ) WindowSizeI *= StartN;    /* geometric */
  288.     else           WindowSizeI += DeltaN;    /* linear */
  289.     }
  290.  
  291.     /* --- compute regression line --- */
  292.  
  293.     RegressI( &SlopeD, &InterceptD, lognFP, logrsFP, CounterI );
  294.     hurstFP[0] = pow( 2., (2.*SlopeD-1.) ) - 1.;    /* full regression */
  295.     for( wxI = 1; wxI < CounterI-(HREG-2); wxI++ ) {
  296.     RegressI( &SlopeD, &InterceptD, &lognFP[wxI-1], &logrsFP[wxI-1], HREG );
  297.     hurstFP[wxI] = pow( 2., (2.*SlopeD-1.) ) - 1.;    /* 3-point regression */
  298.     }
  299.  
  300. Done:
  301.     if ( wxFP  != (float *)0 ) free( (void *)wxFP  );
  302.     if ( rswFP != (float *)0 ) free( (void *)rswFP );
  303.     return( rcI );
  304. }
  305.  
  306. /************************************************************************
  307.  *                                    *
  308.  *    CorrDimI() - Correlation Dimension                *
  309.  *                                    *
  310.  ************************************************************************
  311.     NOTE: Each of the return vectors MUST be of length NIterI or greater.
  312.  
  313.     For each iteration, using the embedding dimension DimI, and
  314.     re-constructing the phase space with a lag of tauI, the
  315.     distance, correlation dimension, log(distance) and log(cd) are
  316.     computed and stored.
  317.  
  318.     The regression coefficients for log(distance)[x] & log(cd)[y] is
  319.     computed for all of the points in the iteration as well as 
  320.     groups of 3 points at a time.
  321.  
  322.     NOTE: the accompanying wingz scripts comptue the correlation
  323.     dimension for a range of embedding dimensions.  This routine
  324.     computes it for only a single embedding dimension.
  325.  */
  326.  
  327. int CorrDimI( rFP, crFP, logrFP, logcrFP, regFP,
  328.           xFP, nI, DimI, tauI, dtD, RD, NIterI )
  329. float    *rFP;        /* (r) distance */
  330. float    *crFP;        /* (r) correlation dimension */
  331. float    *logrFP;    /* (r) log (distance) */
  332. float    *logcrFP;    /* (r) log (correlation dimension) */
  333. float    *regFP;        /* (r) regression of x=log(dist) vs y=log(cd) */
  334. register float    *xFP;    /* (i) input vector */
  335. int     nI;        /* (i) length of input vector */
  336. int     DimI;        /* (i) embedding dimension */
  337. int     tauI;        /* (i) time lag for re-constructing phase space */
  338. double     dtD;        /* (i) change in distance each iteration */
  339. double     RD;        /* (i) initial distance */
  340. int     NIterI;    /* (i) number of iterations to make */
  341. {
  342.     int         i,j,k;        /* work indicies */
  343.     int         itI, ktI;    /* derived work indicies */
  344.     int         nptI;        /* # of points */
  345.     int         StopI;        /* flag to stop iterating */
  346.     double     D;        /* distance */
  347.     double     RRD;        /* current distance squared */
  348.     double     CRD;        /* correlation dimension */
  349.     long     ThetaL;    /* count of "hits" */
  350.     int         rcI=0;        /* return code */
  351.     int         wxI;        /* work index */
  352.     double     vD;        /* work value */
  353.  
  354.     /* NOTE: the following computations are based on Equation (12.2)
  355.      * in "Chos & Order in the Capital Markets"..
  356.      */
  357.  
  358.     nptI = nI - DimI * tauI;    /* max points in embedding dimension */
  359.     for( StopI = 0; StopI < NIterI; StopI++ ) {
  360.     ThetaL = 0;        /* # of phase-space points in radius RD */
  361.     RRD = RD * RD;        /* square of current distance */
  362.     /* NOTE: use this rather than computing the square-root:
  363.      *  if a > b  ->  a*a > b*b   where a>0 && b>0.
  364.      */
  365.     for( k = 0; k < nptI; k++ ) {
  366.         /* NOTE: take advantage of the fact that we are computing
  367.          * the square of the difference between ZFP[j,k] & ZFP[j,i].
  368.          * We need only compute the upper part of a triangular matrix
  369.          * and multiply ThetaL by 2 at the end.
  370.          */
  371.         for( i = k+1; i < nptI; i++ ) {
  372.         D = 0.;                /* distance */
  373.         itI = i;            /* work surrogate for i */
  374.         ktI = k;            /* work surrogate for k */
  375.         for( j = 0; j < DimI; j++ ) {
  376.             /* ZFP[j,i] = xFP[j*tauI + i] */
  377.             /* ZFP[j,k] = xFP[j*tauI + k] */
  378.             /* vD = ZFP[j,k] - ZFP[j,i]; */
  379.             vD = xFP[ktI] - xFP[itI];    /* distance */
  380.             D += vD * vD;        /* squared distance */
  381.             itI += tauI;        /* update phase dimension */
  382.             ktI += tauI;
  383.         }
  384.         if ( D < RRD ) ThetaL++;    /* within the radius */
  385.         }
  386.     }
  387.     vD = ((double)nptI-1.) * ((double)nptI); /* adjust for i <> k */
  388.     CRD = 2.*((double)ThetaL) / vD;    /* correlation dimension */
  389.     /* NOTE: factor of 2 to compensate for only computing upper
  390.      * triangular part of the matrix
  391.      */
  392.     rFP[StopI]    = RD;        /* current distance */
  393.     crFP[StopI]    = CRD;        /* current correlation dimension */
  394.     logrFP[StopI]    = RD>EMIN?  log(RD):0.;        /* log(distance) */
  395.     logcrFP[StopI]    = CRD>EMIN? log(CRD):0.;    /* log(cd) */
  396.  
  397.     RD += dtD;            /* increment the distance */
  398.     }
  399.  
  400.     /* --- Compute the slope --- */
  401.  
  402.     RegressI( &D, &vD, &logrFP[0], &logcrFP[0], StopI );
  403.     regFP[0] = D;
  404.     for( StopI = 1; StopI < NIterI-1; StopI++ ) {
  405.     RegressI( &D, &vD, &logrFP[StopI-1], &logcrFP[StopI-1], 3 );
  406.     regFP[StopI] = D;
  407.     }
  408.     regFP[StopI] = 0.;            /* make things clean */
  409.  
  410. Done:
  411.     return( rcI );
  412. }
  413.  
  414. /************************************************************************
  415.  *                                    *
  416.  *    CorrDim2I() - Alternate Computation of Correlation Dimension    *
  417.  *                                    *
  418.  ************************************************************************
  419.     NOTE: Each of the return vectors MUST be of length NIterI or greater.
  420.  
  421.     For each iteration, using the embedding dimension DimI, and
  422.     re-constructing the phase space with a lag of tauI, the
  423.     distance, correlation dimension, log(distance) and log(cd) are
  424.     computed and stored.
  425.  
  426.     The regression coefficients for log(distance)[x] & log(cd)[y] is
  427.     computed for all of the points in the iteration as well as 
  428.     groups of 3 points at a time.
  429.  
  430.     This program automatically computes the minimum radius (RD) and
  431.     the delta radius (dtD) based on the minimum radius, maximum radius,
  432.     and the number of iterations.
  433.  */
  434.  
  435. #define    RCORR    5    /* regression correlation */
  436.  
  437. int CorrDim2I( rFP, crFP, logrFP, logcrFP, regFP,
  438.         xFP, nI,
  439.         DimI, tauI, NIterI )
  440. float    *rFP;        /* (r) distance */
  441. float    *crFP;        /* (r) correlation dimension */
  442. float    *logrFP;    /* (r) log (distance) */
  443. float    *logcrFP;    /* (r) log (correlation dimension) */
  444. float    *regFP;        /* (r) regression of x=log(dist) vs y=log(cd) */
  445. register float    *xFP;    /* (i) input vector */
  446. int     nI;        /* (i) length of input vector */
  447. int     DimI;        /* (i) embedding dimension */
  448. int     tauI;        /* (i) time lag for re-constructing phase space */
  449. int     NIterI;    /* (i) number of iterations to make */
  450. {
  451.     int         i,j,k;        /* work indicies */
  452.     int         itI, ktI;    /* derived work indicies */
  453.     int         nptI;        /* # of points */
  454.     double     D;        /* distance */
  455.     double     RD;        /* current raidus */
  456.     double     dtD;        /* increment in radius */
  457.     double     CRD;        /* correlation dimension */
  458.     long     ThetaL, wL;    /* count of "hits" */
  459.     int         rcI=0;        /* return code */
  460.     int         wxI;        /* work index */
  461.     double     vD;        /* work value */
  462.     int         FirstFlagI;    /* first time through flag */
  463.     double     RMinD;        /* minimum distance */
  464.     double     RMaxD;        /* maximum distance */
  465.     double     RSclD, ROffD;    /* scale & offset for sorting */
  466.     long    *ThetaLP=0;    /* work array */
  467.     int         WIterI;    /* work iterations */
  468.     int         WMinXI;    /* minimum index */
  469.  
  470.     if ( NIterI < 10 ) { return( -1 ); }
  471.     nptI = nI - DimI * tauI;    /* max points in embedding dimension */
  472.     ThetaLP = (long *)malloc( NIterI*sizeof(long) );
  473.     if ( ThetaLP == (long *)0 ) {
  474.     rcI = -1;
  475.     goto Done;
  476.     }
  477.     memset( (void *)ThetaLP, 0, NIterI*sizeof(long) );
  478.     for( wxI = 0; wxI < NIterI; wxI++ ) {
  479.     rFP[wxI]    = 0.0;
  480.     crFP[wxI]    = 0.0;
  481.     logrFP[wxI]    = 0.0;
  482.     logcrFP[wxI]    = 0.0;
  483.     regFP[wxI]    = 0.0;
  484.     }
  485.  
  486.     /* --- Compute the Min/Max Distance between points --- */
  487.  
  488.     RMinD = RMaxD = 0;
  489.     FirstFlagI = 1;
  490.     for( k = 0; k < nptI; k++ ) {
  491.     for( i = k+1; i < nptI; i++ ) {
  492.         D = 0.;            /* distance to point in phase space */
  493.         itI = i;            /* work surrogate for i */
  494.         ktI = k;            /* work surrogate for k */
  495.         for( j = 0; j < DimI; j++ ) {    /* phase space on the fly */
  496.         vD = xFP[ktI] - xFP[itI];    /* distance */
  497.         D += vD * vD;            /* squared distance */
  498.         itI += tauI;            /* update phase dimension */
  499.         ktI += tauI;
  500.         }
  501.         D = sqrt( D );    /* euclidean distance */
  502.         if ( FirstFlagI ) {        RMinD = RMaxD = D; FirstFlagI = 0; }
  503.         else if ( D < RMinD )    RMinD = D;
  504.         else if ( D > RMaxD )    RMaxD = D;
  505.     }
  506.     }
  507.  
  508.     /* --- Compute the Scale / Offset to sort into bins --- */
  509.  
  510.     if ( (RMaxD - RMinD) < EMIN ) {
  511.     rcI = -2;
  512.     goto Done;
  513.     }
  514.     WIterI = 3 * NIterI;    /* part of space to explore */
  515.     WMinXI = NIterI;        /* first part to ignore */
  516.     RSclD = ((double)WIterI) / (RMaxD - RMinD);
  517.     ROffD = -RSclD * RMinD - WMinXI;
  518.  
  519.     /* --- Sort everything into bins --- */
  520.  
  521.     for( k = 0; k < nptI; k++ ) {
  522.     for( i = k+1; i < nptI; i++ ) {
  523.         D = 0.;            /* distance to point in phase space */
  524.         itI = i;            /* work surrogate for i */
  525.         ktI = k;            /* work surrogate for k */
  526.         for( j = 0; j < DimI; j++ ) {    /* phase space on the fly */
  527.         vD = xFP[ktI] - xFP[itI];    /* distance */
  528.         D += vD * vD;            /* squared distance */
  529.         itI += tauI;            /* update phase dimension */
  530.         ktI += tauI;
  531.         }
  532.         D = sqrt( D );    /* euclidean distance */
  533.         wxI = RSclD * D + ROffD;            /* index into table */
  534.         if ( wxI <  0      ) wxI = 0;
  535.         if ( wxI <  NIterI ) ThetaLP[wxI] += 2;    /* acct for full mtx */
  536.     }
  537.     }
  538.  
  539.     /* --- Compute the Correlation Dimension --- */
  540.  
  541.     vD = ((double)nptI - 1.) * ((double)nptI);    /* adjust for i <> k */
  542.     ThetaL = 0;
  543.     dtD = (RMaxD - RMinD) / ((double)WIterI);    /* space between bins */
  544.     RD  = RMinD + (WMinXI+1) * dtD;        /* first radius */
  545.     for( wxI = 0; wxI < NIterI; wxI++ ) {
  546.     ThetaL += ThetaLP[wxI];        /* cumulative count */
  547.     CRD = ((double)ThetaL) / vD;    /* correlation dimension */
  548.     rFP[wxI]    = RD;        /* current distance */
  549.     crFP[wxI]    = CRD;        /* current correlation dimension */
  550.     logrFP[wxI]    = RD>EMIN?  log(RD):0.;        /* log(distance) */
  551.     logcrFP[wxI]    = CRD>EMIN? log(CRD):0.;    /* log(cd) */
  552.     RD += dtD;            /* increment the distance */
  553.     }
  554.  
  555.     /* --- Compute the slope --- */
  556.  
  557.     RegressI( &D, &vD, &logrFP[0], &logcrFP[0], NIterI );
  558.     regFP[0] = D;
  559.     for( wxI = 1; wxI < NIterI-(RCORR-2); wxI++ ) {
  560.     RegressI( &D, &vD, &logrFP[wxI-1], &logcrFP[wxI-1], RCORR );
  561.     regFP[wxI] = D;
  562.     }
  563.     regFP[wxI] = 0.;            /* make things clean */
  564.  
  565. Done:
  566.     if ( ThetaLP != (long *)0 ) free( (void *)ThetaLP );
  567.     return( rcI );
  568. }
  569.  
  570. /************************************************************************
  571.  *                                    *
  572.  *    LyapunovI() - Compute the Largest Lyapunov Exponent        *
  573.  *                                    *
  574.  ************************************************************************
  575.     The return arrays are:
  576.     zlyapFP        - lyapunov exponent estimate
  577.     EvolveItFP    - Evolution * Iterations
  578.     diFP        - distance
  579.     DFFP        - starting distance
  580.     regFP        - regression [0=all; 1.. = 3-element]
  581.     Each of these MUST be dimensioned nI/EvolveI in length.
  582.  
  583.     The function returns negative return codes for errors.
  584.     Positive return codes are the number of items in the return
  585.     arrays.
  586.  
  587.     This code has been verified against Wolfe's FORTRAN program.
  588.  */
  589.  
  590. int LyapunovI(  zlyapFP, EvolveItFP, diFP, DFFP,
  591.         xFP, nI, DimI, tauI, dtD, ScaleMaxD, ScaleMinD,
  592.         EvolveI, LagI )
  593. float    *zlyapFP;    /* (r) lyapunov exponent estimate */
  594. float    *EvolveItFP;    /* (r) evolution * iterations */
  595. float    *diFP;        /* (r) distance at this level */
  596. float    *DFFP;        /* (r) maximum distance at start */
  597. float    *xFP;        /* (i) input data vector */
  598. int     nI;        /* (i) number of points in data vector */
  599. int     DimI;        /* (i) embedding dimension */
  600. int     tauI;        /* (i) phase shift for reconstructing phase space */
  601. double     dtD;        /* (i) increase in each measurement */
  602. double     ScaleMaxD;    /* (i) maximum divergence */
  603. double     ScaleMinD;    /* (i) minimum divergence */
  604. int     EvolveI;    /* (i) number of iterations to evolve */
  605. int     LagI;        /* (i) minimum time between pairs */
  606.             /*     NOTE: tauI <= LagI < 12. */
  607. {
  608.     int         nptI;        /* effective number of points in data set */
  609.     int         ind1I, ind2I;    /* indices into phase */
  610.     int         indoldI;    /* last ind2I value */
  611.     int         wI;        /* work integer */
  612.     int         i,j,k;        /* work indices for loops */
  613.     int         jTauI;        /* work variable to reduce computation */
  614.     int         FirstFlagI;    /* first time flag */
  615.     int         rcI=0;        /* return code */
  616.     double     D;        /* work value */
  617.     double     diD;        /* "di" in wingz */
  618.     double     diiD;        /* "dii" in wingz */
  619.     double     vD;        /* temp work value */
  620.     double     DF;        /* final divergence */
  621.     double     TotalD;    /* accumulator */
  622.     int         ItersI;    /* # of iterations */
  623.     double     zlyapD;    /* lyapunov exponent? */
  624.     double     zmultD;    /* z-multiplier */
  625.     double     anglmxD;    /* maximum angle */
  626.     double     thminD;    /* theta minimum */
  627.     double     DNewD;        /* new distance */
  628.     double     dotD;        /* ? */
  629.     double     cthD;        /* cosine theta */
  630.     double     thD;        /* theta */
  631.     double    *pt1DP=0;    /* pointer to work area */
  632.     double    *pt2DP=0;    /* ditto... */
  633.     double     Log2D;        /* log of 2.0 */
  634.     double     mD, bD;    /* slope & intercept */
  635.  
  636.     if ( (pt1DP=(double *)malloc( (DimI+2)*sizeof(*pt1DP) )) == (double *)0
  637.       || (pt2DP=(double *)malloc( (DimI+2)*sizeof(*pt2DP) )) == (double *)0 ) {
  638.     rcI = -1;
  639.     goto Done;
  640.     }
  641.  
  642.     for( wI = 0; wI < nI; wI++ ) {
  643.     zlyapFP[wI]    = 0.0;
  644.     EvolveItFP[wI]    = 0.0;
  645.     diFP[wI]    = 0.0;
  646.     DFFP[wI]    = 0.0;
  647.     }
  648.  
  649.     nptI = nI - (DimI * tauI + EvolveI);
  650.  
  651.     FirstFlagI = 1;            /* first time hit */
  652.     ind1I = ind2I = 0;            /* first index */
  653.     diiD = diD = 0.;
  654.     for( i = LagI+1; i < nptI; i++ ) {
  655.     D = 0.;
  656.     for( j = jTauI = 0; j < DimI; j++, jTauI += tauI ) {
  657.         vD  = xFP[ jTauI + ind1I ] - xFP[ jTauI + i ];
  658.         D  += vD * vD;
  659.     }
  660.     D = sqrt(D);
  661.     if ( (FirstFlagI || D <= diD) && D >= ScaleMinD ) {
  662.         diD = D;        /* smallest point > ScaleMinD */
  663.         ind2I = i;
  664.         FirstFlagI = 0;
  665.     }
  666.     }
  667.     if ( FirstFlagI ) { rcI = -2; goto Done; }
  668.  
  669.     TotalD = 0.0;
  670.     Log2D = log(2.);
  671.     for( ItersI = 1; ; ItersI++) {
  672.     /* --- Save Coordinates of Evolved points --- */
  673.     DF = 0.;
  674.     for( j = 0; j < DimI; j++ ) {
  675.         pt1DP[j] = xFP[ ind1I + EvolveI + j * tauI ];
  676.         pt2DP[j] = xFP[ ind2I + EvolveI + j * tauI ];
  677.         vD = pt1DP[j] - pt2DP[j];
  678.         DF += vD * vD;
  679.     }
  680.     DF = sqrt( DF );    /* final divergence */
  681.     TotalD += log(DF/diD) / (EvolveI * dtD * Log2D );
  682.     zlyapD = TotalD / ItersI;
  683.  
  684.     /* --- pass data back to calling program --- */
  685.  
  686.     wI = ItersI - 1;
  687.     zlyapFP[wI]    = zlyapD;
  688.     EvolveItFP[wI]    = ((double)EvolveI) * ((double) ItersI);
  689.     diFP[wI]    = diD;
  690.     DFFP[wI]    = DF;
  691.     if ( diD > 5000. || DF > 5000. )
  692.         break;
  693.  
  694.     indoldI = ind2I;
  695.     zmultD = 1.0;
  696.     anglmxD = 0.3;
  697.     for(;;) {
  698.         thminD = 3.14;
  699.         /* --- look for replacement points --- */
  700.         for( i = 0; i < nptI; i++ ) {
  701.         wI = i - (ind1I + EvolveI);
  702.         if ( -LagI < wI && wI < LagI )
  703.             continue;        /* too close */
  704.         DNewD = 0.;
  705.         for( jTauI=i, j=0; j < DimI; j++, jTauI += tauI ) {
  706.             vD = pt1DP[j] - xFP[jTauI];
  707.             DNewD += vD * vD;
  708.         }
  709.         DNewD = sqrt(DNewD);
  710.         if ( DNewD > (zmultD * ScaleMaxD) || DNewD < ScaleMinD )
  711.             continue;
  712.  
  713.         dotD = 0.;
  714.         for( jTauI=i, j=0; j < DimI; j++, jTauI += tauI )
  715.            dotD += (pt1DP[j] - xFP[jTauI]) * (pt1DP[j] - pt2DP[j]);
  716.         cthD = fabs( dotD / (DNewD * DF) );    /* per Peters */
  717.         if ( cthD > 1. ) cthD = 1.;
  718.         thD = acos(thD);    /* thD = cos(thd) does not make sense */
  719.         if ( thD <= thminD ) {
  720.             thminD = thD;
  721.             ind2I  = i;
  722.             diiD   = DNewD;
  723.         }
  724.         }
  725.  
  726.         if ( thminD < anglmxD ) goto EvolveInd1;
  727.         /* --- thminD >= anglmxD --- */
  728.         zmultD += 1.;
  729.         if ( zmultD < 5. ) continue;
  730.         zmultD = 1.;
  731.         anglmxD *= 2.;
  732.         if ( anglmxD >= 3.14 )
  733.         break;
  734.     }
  735.     ind2I = indoldI + EvolveI;
  736.     diiD = DF;
  737.     EvolveInd1:
  738.     ind1I += EvolveI;
  739.     if ( ind1I >= nptI ) break;
  740.     diD = diiD;
  741.     }
  742.     rcI = ItersI;
  743.  
  744. Done:
  745.     if ( pt1DP != (double *)0 ) free( (void *)pt1DP );
  746.     if ( pt2DP != (double *)0 ) free( (void *)pt2DP );
  747.     return( rcI );
  748. }
  749.